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Abstract 

We re-analyse the data on J/ip cross sections for proton-nucleus 
and nucleus-nucleus collisions obtained by the CERN experiments 
NA38, NA51 and NA50. Our systematic analysis of error propagation 
shows that the no anomalous suppression of J/ip in Pb-Pb collisions 
can be substantiated at 95% confidence limit. The significance of the 
NA50 result would be enhanced by more precise p-A data on J/ip 
cross sections. 
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Recently the NA50 experiment at CERN announced an anomalous sup- 
pression of J /if) in Pb-Pb collisions at a beam energy of 158 GeV/nucleon 
[|I], 0. Such an observation had earlier been proposed as a signal for the 
formation of a quark-gluon plasma (QGP) 0. In view of the implications 
of this observation, we reanalyse the CERN data in this paper, taking into 
account many obvious sources of uncertainties. We conclude that the data 
do not yet bear the claim they are asked to support. 

The NA50 experiment as well as the NA51 and NA38 experiments have 
taken data on J /if) production in proton-nucleus collisions with beam energies 
of 450 and 200 GeV, and in nucleus-nucleus collisions at beam energies of 
200 and 158 GeV/nucleon. The CM energies in equivalent proton-proton 
collisions are y/S = 29.1, 19.4 and 17.2 GeV. The acceptance region, 5, for 
the data analysed in is < y < 1 (where y is the rapidity of the J / if) 
in the CM frame of colliding protons at the same v^S*) and | cos6*\ < 0.5 
(where 8* is the angle of the decay muon in the Collins-Soper frame). 

The analysis consists of the following steps — 

1. All the data on the cross section per nucleon (a/AA') are rescaled 
to 200 GeV beam energy, using a fit to previous data on the energy 
dependence of this cross section. The 158 Gev data are rescaled by a 
factor 1.32 and the 450 GeV data by the factor 0.48. 

2. The rescaled data, at a common v^S* = 19.4 GeV, are extrapolated to 
| cos8*\ < 1, through a multiplicative factor of 2. 

3. The resulting cross section is extrapolated to the full forward hemi- 
sphere, y > 0, or equivalently, x F > by multiplying by 1.07. 

4. The rescaled pA data in the forward hemisphere are fitted to a power- 
law 

Ba(A) = a A a , (1) 

where B is the branching ratio for J /if) — > and A is the effec- 

tive mass number, given by the product of the mass numbers of the 
projectile and target. 

It is found that the data obtained from light ions (O-Cu, O-U and S-U) 
lie on the fitted line, whereas the Pb-Pb data lie far below. While the full 
procedure outlined above is used in 0, the steps 2 and 3 are not performed 
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in 0. However, the extrapolation in y/S in both [M 0] uses a fit to the world 
data for x F > and full coverage in cos 9*. 

We note that there are uncertainties in every step of this procedure, and 
also in the order in which these steps are to be performed. Each such un- 
certainty is a source of error which must be combined with the statistical 
and systematic errors in the Pb-Pb data. We present three seperate analy- 
ses. First we accept the extrapolation in step 1 of the NA50 analysis |jj and 
examine only the error in step 4 above. The Pb-Pb data is seen to be con- 
sistent with the p-A data at the 95% confidence level (CL). Next we analyse 
the errors due to steps 1 and 4 only. Note that when the extrapolations are 
performed by NA50 |], the errors in the extrapolations in steps 2 and 3 are 
not needed in estimating the confidence limits of the final result. Again the 
p-A and Pb-Pb data are consistent at the 95% CL. Finally, we examine the 
situation when steps 2 and 3 are performed before step 1, and combine the 
errors in each of these four steps. The p-A and Pb-Pb data are once again 
seen to be consistent at the 95% CL. 

Since we shall repeatedly have to propagate errors, we outline the usual 
statistical procedure for this ||. If the expectation values of a set of variables, 
Pi, are known along with their full covariance matrix, c^, then the errors in 
any function of these variables, f{pi) can be easily written down. Recall 
that the diagonal elements of the covariance matrix are simply related to the 
errors, Api = ^/c^. The off-diagonal elements can be used to extract the 
normalised correlation coefficients = —c^j Jc^cfj. These lie in the range 
\tij\ < 1. If f(pi) is a linear function of the pi, then (f(pi)) = f((Pi)), and 
the error 

(a/) 2 = Ee## ( 2 ) 

Note that when the correlations vanish, this reduces to the usual formula for 
adding errors in quadrature^. In the succeeding analysis it will turn out that 
we often have information on the full covariance matrix. As a result, we can 
easily evaluate the formula in eq. @. In the exceptional cases when we have 
no knowledge of the off-diagonal elements, we will choose them (by choosing 



3 When / is not a linear function, the expectations and errors can be obtained either 
under the assumption of a multivariate Gaussian distribution of the pi, or by a detailed 
bootstrap or Monte Carlo Q. In this paper we will not need to use a non-linear function 
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Tij) to minimise the errors. This would give the best possible statistical 
significance to the result. 

We begin by examining the errors involved in step 4 above. The form in 
eq. (|l]) is linearised by taking the logarithm of both sides. A fit is performed 
using the p-A data listed in 0] obtained by rescaling to a common VS = 19.4 
GeV. The result of this fit is — 

a = 2.33(1±0.11) nb, a = 0.908±0.026, Cov(log(x , a) = -0.0026. (3) 

Note that the value of a is exactly the same as quoted [[[], 0] by NA50. 
Following exactly the NA50 analysis in Figure |l] we show the best fit 
line to the pA data along with the error bands about this fit line (using eq. 
|). The prediction for Pb-Pb collisions is Ba(Pb-Pb) = 0.87(1 ± 0.18) nb. 
The measured point, Ba(Pb-Pb) = 0.59 ± 0.04 nb, is then within 2a of the 
extrapolation. Note that the extrapolation errors are more crucial than the 
statistical errors on the Pb-Pb cross section. As a result, the J/ip suppression 
observed by NA50 is indicative, but not a quantitative proof of the existence 
of new physics, such as the QGP. 

The assumption that O-Cu, O-U and S-U collisions do not yield a QGP 
can be used to make a similiar statement. With this assumption, these three 
extra points are included in the fit. Then we find 

(x = 2.31(1 ±0.08) nb, a = 0.911 ±0.014, Cov(loga , a) = -0.001. (4) 

The prediction for Pb-Pb collisions is Ba(Pb-Pb) = 0.89(1 ±0.09) nb. The 
measured point is then less than 3a away from the prediction. The extra 
assumption makes the observation stronger. 

We turn now to the second analysis. Recall that the rescaling of energy is 
performed using a fit to data on the Vo dependence of cross sections in the 
full forward hemisphere, using a formula discussed later. We have no reason 
to believe that the two ratios 



Ba(y> 0, |cos0*| <1,VS) 



r 



= Ba(0<y<l,\cos6*\<l,VS) 
' ' N ' Ba(0 <y <l,\cos6*\<l,Vs') 



(5) 



are equal. The higher the energy, the larger is the allowed range of y, and the 
slice < y < 1 is consequently a smaller part of the full phase space. For the 
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Figure 1: Data on the J ftp cross section for < y < 1 and | cos6>*| < 0.5, 
plotted against the product of the mass numbers of the target and projectile. 
The best fit power-law for the nuclear effects is shown along with the la 
(enclosed by full lines), and the 2a (dashed line). The error bar on the 
Pb-Pb point shows the 1 and 2a errors with the same line styles. 



same reason, the orders in which the two extrapolations are performed may 
give different results. However, cuts on cos6** are not expected to change 
these ratios provided they are identical for the numerator and denominator. 

Since we do not have data on the dependence of Bas, we use models 
to extract r. We have at our disposal the colour evaporation model || and 
the colour octet model @. Both give a reasonably good description of the 
\/S dependence of Ba y>0 as well as the distribution in y (or x F ). For details 
we refer the readers to || |9|. In Figure we show the ratio r obtained in 
these models using several different sets of modern structure functions. Note 
the good agreement among the models and the structure function sets. The 
results shown in the figure translate to 

r(19.4,17.2) = 1.382±°;°it and r(19.4,29.1) = 0.412™, (6) 

where the central values are the averages over the model predictions and the 
errors are the la range of errors under variation of the structure function 
set. These ratios and their errors may also be obtained from data if they are 
available. 

If the procedure followed by the NA50 group is adopted, then the errors 
from the recaling of energy must be compounded with the statistical errors 
on the data. We can add these errors in quadrature since the errors in r 
are independent of the statistical errors in the NA50 and NA38 experiments. 
In addition, the errors in the fits to the nuclear effects must be included 
as before. The extrapolation to the full forward phase space is a common 
multiplicative factor for all the data, and is unnecessary for the comparison. 
Again our conclusions will therefore apply to both Refs. and @. Using 
this procedure, we find that the extrapolation of the p-A data to Pb-Pb is 
1.0(1±0.21) nb, whereas the measured value after rescalings is 0.61 ± 0.04 
nb. Thus the Pb-Pb data is within 2cr of the extrapolation from the p-A 
data. 

We move on to the third analysis. Recall that the Collins-Soper frame is 
the rest frame of the decaying particle with the z-direction chosen such that 
the spin is aligned in this direction if the helicity is positive. Then, in this 
frame, the angular distribution of the decay leptons coming from any massive 
vector particle (either 7* for the Drell-Yan process or any vector quarkonium) 
should be 1 + cos 2 8*, when strong interaction corrections are neglected. Then 
the extrapolation from | cos#*| < 0.5 to the full acceptance is a factor 32/13 
for both Drell-Yan and vector quarkonium production, independent of y/S. 
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Figure 2: The scaling factor r (eq. |5|) in the colour octet model (full circles) 
and the colour evaporation model (open circles) for various structure func- 
tions. The upper set is for r(19.4, 17.2), and the lower for r(19.4, 29.1). The 
error bars show the variation due to varying all model parameters within 
one standard deviation and the charm quark mass between 1.4 GeV and 1.6 
GeV. 
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The NA50 group uses this factor for Drell-Yan cross sections, but a factor 
2 for J/ip. Since this is an overall factor, it does not affect the present anal- 
ysis. However, in any comparison of Drell-Yan and J /if) cross sections, the 
correct factorfj must be used. It should also be recalled that QCD corrections 
change this distribution. The modification can be written as a series in the 
strong-coupling, a s , multiplying cos 2 #*. The first term in this series is 1 for 
all vector particles, but the subsequent terms may be different for 7* and 
J/ip. Then the extrapolation uncertainty in the ratio of Drell-Yan and J /if) 
cross sections is of the order of a s . 

There remains the question of extrapolation from 0<y<ltoy>0 
(or equivalently, x F > 0). Many experiments [n| have parametrised the x P 
distribution in the form (1 — x F ) d . Using this form, and the range of errors 
quoted, we find that 





1.03- 


1.05 


(Vs = 


16.8) 


Ba(y > 0) = Ba(0 < y < 1) x < 


1.06- 


1.07 


(Vs = 


20.4) 




1.11 - 


1.26 


(y/s = 


27.4) 



(7) 



In order to obtain these scalings at the values of \/S we need, we could use a 
parametrisation of the \/~S dependence of the exponent d fitted in the form 







d 



(8) 



1 + b/y/S' 

However, this gives results at variance with those in eq. (|7|). This might be 
due to fitting too simple a functional form to the data. On the other hand, 
the colour octet model calculations yield numbers in reasonable agreement 
with eq. ([?]) and the actual value used by the NA50 analysis at y/S = 19.4 

4 After this paper was written we were informed by some members of the NA50 collab- 
oration that data on J/ip production at lower energies shows a flat angular distribution, 
suggesting 2 as the appropriate rescaling factor. 
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GeV. We find 



Ba(y > 0) = Ba(0 < y < 1) x 



1.053 ±0.009 (VS = 17.2) 

1.082 ±0.014 (^=19.4) (9) 

1.224 ±0.041 (^ = 29.1) 



The third analysis may be performed by doing the extrapolation from 
5 to the forward hemisphere first, and then using the available data on the 
energy dependence of Ba y> o to bring all the data to a common value of y/S. 
In this case the errors in eq. (§) must be combined with the error in energy 
extrapolation. 

The rescaling of Ba with v^S* is performed by a fit of the form (Tl[ 



Ba(VS) = o-o f 1 



M 



j/i> 



S 



(10) 



to data on the total forward J / ip cross section at different v^- The problem 
is linearised by taking logarithms on both sides of this equation. Taking 
Mi 



i W = 3.097 GeV, JTTJ reports the best fit to bd3 
38(1 ±0.16) nb, 



<7o 



and 



n 



12.0 ±0.9. 



(11) 



These numbers were obtained by fitting both p-p and p-A data. Nuclear 
effects in the p-A data were eliminated by assuming a = 0.91 in The 
formula in eq. (|l~0l), with the values in eq. (|TT|), was used by the NA50 group 
to compute the rescaling in y/S fl], [2|. We find 

#(19.4, 17.2) = 1.34 ±0.03 and #(19.4,29.1) = 0.48 ±0.03, (12) 

where the error in R depends only on the error in n. Possible systematic 
errors in this extrapolation can be investigated by changing the form in eq. 
(p!U|). In particular, at larger energies, one should take into account the 
expected rise in cross sections as (roughly) log S. 

5 Using the data collected in Table 8 of |llj and using MINUIT without linearising 
the problem, we get somewhat different values of the parameters — tro = 44^}{ n b and 



n 



13.0t^. This evidently increases the errors quoted below on the quantity R. 
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Figure 3: Data on the J/ip cross section Ba(y > 0, 19.4) plotted against the 
product of the mass numbers of the target and projectile. The best fit power- 
law for the nuclear effects is shown along with the 1 and 2a error bands (full 
and dashed lines respectively). Similiarly, the 1 and 2a errors on the Pb-Pb 
point are also shown. 
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Using eqs. (|9|) and (0), along with a scaling factor of 32/13 for the 
extrapolation of the Collins-Soper angle, we find that 



Ba(y>0,\co89*\ < 1,19.4) = (1.41 ± 0.02) x Ba 5 (l7.2) 
Ba(y > 0, | cos 6* | < 1, 19.4) = (0.59 ± 0.02) x Ba s (29.1). 

Recall that the factor of 32/13 is a common overall factor, and hence its value 
or error is immaterial to the following analysis. The errors in the coefficients 
in eq. (p~3| ) are independent of the measurements of the NA38, NA50 and 
NA51 experiments. Hence they may be added in quadrature to the statistical 
errors on the data. When the pA data are scaled in this fashion, the fit of 
the form in eq. ([T]) becomes 

cr = 7.45(1 ±0.11) nb, a = 0.882 ± 0.026, Cov(log<r , a) = -0.0028. 

(14) 

Although the value of a decreases, it still remains consistent with the E-789 
extraction of a = 0.90 ± 0.02 [Q. The scaling of the CERN data, however, 



is a little stronger than that observed by the E-772 experiment^ which sees 



a = 0.920 ±0.008 fll3[ . The comparison of the extrapolation of this data with 
the Pb-Pb cross section is shown in Figure |3|. Note that the cross section 
for J /if) production in Pb-Pb collisions is still consistent with those for p-A 
collisions at the 2a level. 

If the light-ion data is included in the fit to the form in eq. ([!]) , by making 
the assumption that there is no new physics seen with light-ions, then the 
fitted parameters become 

a = 7.23(1 ± 0.08) nb, a = 0.890 ± 0.014, Cov(log<7 , a) = -0.0011. 

(15) 

With this extra assumption, the extrapolated line fails to agree with the 
Pb-Pb data at the 2a level, but remains consistent at the 3a level. 

The experiments have more information than just the total cross section. 
In particular, they measure the transverse energy, E T , dependence of the 
J/i/j cross section. A quantitative analysis Jl4| found that the well-known 



absorption mechanism can account for the E T dependence in S-U collisions, 



6 The comparison of all these values of a is a measure of the systematic uncertainty in 
this quantity. One can assign the value 0.02 to the systematic error on a. 



10 



but not in Pb-Pb collisions. This analysis uses the survival probability of the 
J/if), defined by the ratio 



a AA '(y>0 7 


cos 0*| <l,y/S= 19.4, E T ) 


aP D (y > 0, 


cos 8* 


< 1, y/S = 19.4, E T ) 



where the cross sections have been scaled to a common y/S and the full 
forward hemisphere, by the prescription of NA50 0. It was found that, 
for Pb-Pb collisions, the data point at the lowest E T agrees with a Glauber 
description, but the remaining points fall below the prediction [TJJ . 

This analysis could also be redone in a fashion similiar to our analysis of 
the total forward cross sections. Figure |] shows the data points which are 
scaled as in the second analysis (eq. |6|). All data agree with the Glauber 
model at the la level. The fit to the S-U data is not as good as before; 
although remaining consistent at the la level it seems to lie systematically 
above the Glauber model prediction. The Pb-Pb data agree much better, 
with the sole exception of the point at highest E T . This point agrees with 
the model only at the 2a level. When eq. flTB| ) is used for the rescalings, 
both the S-U and Pb-Pb data lie systematically below the Glauber model 
predictions, although being consistent at the 2a level. In either case, the 
Glauber model does equally well (or badly) in explaining the E T dependence 
of both the S-U and Pb-Pb data. 

In conclusion, we find little statistical significance for an anomalous sup- 
pression of J/ip cross section in Pb-Pb collisions at the CERN SPS energy. 
In our analysis, the uncertainties are dominated by the errors due to the 
extrapolation of the usual nuclear effects in p-A collisions to Pb-Pb. Clearly, 
reduced errors on the p-A cross sections would enhance the statistical signifi- 
cance of the deviations seen by NA50. The importance of controlling nuclear 
effects was first pointed out in |L5| ; it continues to be an important problem. 



Even in future experiments at the RHIC and LHC, the p-A data will always 
serve as a baseline of "known physics" . We feel that the worldwide heavy-ion 
effort must be supplemented by an effective study of the systematics of p-A 
collisions. Unless improved data from p-A collisions are available, interpre- 
tations of data from present and future heavy-ion colliders may be plagued 
by uncertainties of the kind we have pointed out. 
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Figure 4: Data on the E T dependence of the J / ip survival probability S (eq. 
|H3p . The band shows the prediction of a Glauber model [|TJJ within a la 
variation of the model parameters. 
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